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1. Introduction 



This article deals with the following semilinear heat equation with Neu- 
mann boundary conditions: 



u x (l,t) 
u x (Q,t) 
u{x, 0) 



u xx - g\{u) 
ag 2 (u(l,t)) 


uq{x) > 



in (0,1) x [0,T), 
in [0,T), 
in [0,T), 
in [0,1], 



where <7i,<?2 G C 3 (K) are analytic functions in x = and a is a positive 
constant. The nonlinear heat equation models many physical, biological 
and engineering phenomena, such as heat conduction (see, e.g., Ql §20.3], 
@, §1.1]), chemical reactions and combustion (see, e.g., @, §5.5], 4, §1.7]), 
growth and migration of populations (see, e.g., [5|, Chapter 13], 2j, §1.1]), 
etc. In particular, "power-law" nonlinearities have long been of interest as a 
tractable prototype of general polynomial nonlinearities (see, e.g., 0, §5.5], 
0, Chapter 7], 0, B, 0, §1.1]). 

The long-time behavior of the solutions of 
ied (see, e.g., g, 0, Q, 0, ' 



0, Hi] 




has been intensively stud- 
16| and the references 
therein). In order to describe the dynamic behavior of the solutions of ([T]) 
it is usually necessary to analyze the behavior of the corresponding station- 
ary solutions (see, e.g., [l3T ]. 0]), i.e., the positive solutions of the following 
two-point boundary-value problem: 




9l{u) 

a# 2 (n(l)), 
0. 



in (0,1), 



The usual numerical approach to the solution of ([TJ consists of consid- 
ering a second-order finite-difference discretization in the variable x, with 
a uniform mesh, keeping the variable t continuous (see, e.g., 17J]). This 
semi-discretization in space leads to the following initial-value problem: 



u fc (0) 



ui) - gi{ui), 

-£i{u k+ i - 2u k + u k -\) - gi(u k ), (2 < k < n 



uo(x k ), 



u n ) - gi(u n ) + 2 fg2{u n ), 

(1 < k < n) 



(3) 



where h := l/(n— 1) and xi, . . . , x n define a uniform partition of the interval 
[0, 1]. A similar analysis to that in J18f] shows the convergence of the positive 
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solutions of ([3|) to those of ([I]) and proves that every bounded solution of 
([3]) tends to a stationary solution of ([3]), namely to a solution of 



= j?{u2 ~ Ul) - 01 (ui), 

= -^(uk+i - 2u k + u k -i) - gi{v,k), (2<k<n-l) 

= 7^(u n _i - u n ) - g\{u n ) + f^k). 

(4) 

Hence, the dynamic behavior of the positive solutions of (j3J) is rather deter- 
mined by the set of solutions {ui, . . . , u n ) E (R>o) n of (jU). 

Very little is known concerning the study of the stationary solutions of 
(131 and the comparison between the stationary solutions of © and ([1]). In 



[131 ]. [181 ] and [19] there is a complete study of the positive solutions of (jl]) 
for the particular case g\(x) := x p and <?2(#) := x q , i.e., a complete study of 
the positive solutions of 

= Jr(«2 -ui) -u{, 

= jp(u k+ i — 2u k + u k -i) - v%., (2<fc<n-l) (5) 
= ^i{u n -i - u n ) - vPn + ^ful. 

In [13] it is shown that there are spurious solutions of (jlj) for q < p < 2q — 1, 
that is, positive solutions of (jU) not converging to any solution of ([2]) as the 
mesh size h tends to zero. 

In [3] and [3] there is a complete study of (j3|) for p > 2q — 1 and 
p < q. In these articles it is shown that in such cases there exists exactly 
one positive real solution. Furthermore, a numeric algorithm solving a given 
instance of the problem under consideration with operations is pro- 

posed. In particular, the algorithm of [191 ] has linear cost in n, that is, this 
algorithm gives a numerical approximation of the desired solution with 0(n) 
operations. 

We observe that the family of s yste ms has typically an exponential 
number 0(p n ) of complex solutions ([201]). and hence it is ill conditioned from 
the point of view of its solution by the so-called robust universal algorithms 



(cf. [2l(, [22|, |23j])- An example of such algorithms is that of general 
continuation methods (see, e.g., [24]). This shows the need of algorithms 
specifically designed to compute positive solutions of "structured" systems 
like @. 

Continuation methods aimed at approximating the real solutions of non- 
linear systems arising from a discretization of two-point boundary-value 
problems for second-order ordinary differential equations have been con- 
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sidered in the literature (see, e.g., [25(, [26|], [271] )■ These works are usu- 



ally concerned with Dirichlet problems involving an equation of the form 
Uxx = f(x,u,u x ) for which the existence and uniqueness of solutions is 
known. Further, they focus on the existence of a suitable homotopy path 
rather on the cost of the underlying continuation algorithm. As a conse- 
quence, they do not seem to be suitable for the solution of On the 
other hand, it is worth mentioning the analysis of [28] on the complexity of 
shooting methods for two-point boundary value problems. 

Let <7i,<?2 ^ C 3 (IR) be analytic functions in x = such that <?i(0) = 0, 
g[{x) > 0, g'-(x) > and g"'(x) > for all x > with i = 1,2. We observe 
that gi and gi are a wide generalization of the monomial functions of system 
([5]). Moreover, we shall assume throughout the paper that the function 
g := <?i/<?2 is strictly decreasing, generalizing thus the relation p < q in 
([5]). In this article we study the existence and uniqueness of the positive 
solutions of (|4]), and we obtain numerical approximations of these solutions 
using homotopy methods. In a forthcoming paper, we shall consider the 
generalization of the relations q < p < 2q — 1 and 2q — 1 < p in . 

1.1. Our contributions 

In the first part of the article we prove that (|4|) has a unique positive so- 
lution, and we obtain upper and lower bounds for this solution independents 



of h, generalizing the results of 1_9(] 



In the second part of the article we exhibit an algorithm which computes 
an ^-approximation of the positive solution of Such an algorithm is 
a continuation method that tracks the positive real path determined by 
the smooth homotopy obtained by considering Q as a family of systems 
parametrized by a. Its cost is roughly of n log log e arithmetic operations, 
improving thus the exponential cost of general continuation methods. 

The cost estimate of our algorithm is based on an analysis of the condi- 
tion number of the corresponding homotopy path, which might be of inde- 
pendent interest. We prove that such a condition number can be bounded 
by a quantity independent of h := \/n. This in particular implies that each 
member of the family of systems under consideration is significantly better 
conditioned than both an "average" dense system (see, e.g., 0, Chapter 
13, Theorem 1]) and an "average" sparse system ([3fJ, Theorem 1]). 

1.2. Outline of the paper 

Section [2] is devoted to determine the number of positive solutions of 
@. For this purpose, we prove that the homotopy of systems mentioned 
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above is smooth (Theorem [5|). From this result we deduce the existence and 
uniqueness of the positive solutions of (|3J). 

In Section [3] we obtain upper and lower bounds for the coordinates of 
the positive solution of (jU). 

In Section U] we obtain estimates on the condition number of the homo- 
topy path considered in the previous section (Theorem I14p . Such estimates 
are applied in Section [5] in order to estimate the cost of the homotopy con- 
tinuation method for computing the positive solution of (JU). 

2. Existence and uniqueness of stationary solutions 

Let Ui, . . . ,U n be indeterminates over M. Let g\ and g 2 be two functions 
of class C 3 (M) such that gi (0) = 0, g'^x) > 0, g'![x) > and g'{'(x) > for 
all x > with i= 1,2. As stated in the introduction, we are interested in 
the positive solutions of (jU) for a given positive value of a, that is, in the 
positive solutions of the nonlinear system 

' = -{U 2 - U x ) + f ^{Ux), 

< = -(U k+1 -2U k + U k „ 1 ) + h 2 g 1 (U k ), (2<fc<n-l) 
= -(U n ^ 1 -U n )^g 1 (U n )-hag 2 (U n ), 

(6) 

for a given value a = a* > 0, where h := l/(n — 1). Observe that, as a 
runs through all possible values in M>o, one may consider © as a family of 
nonlinear systems parametrized by a, namely, 

' = -{U 3 - Ux) + £<&(E/i), 

< = -{U k+l -2U k + U k ^) + h 2 gi {U k ), (2<k<n-l) 
= -{U n „ l -U n %g 1 {U n )-hAg 2 {U n ), 

(7) 

where A is a new indeterminate. 

2.1. Preliminary analysis 

Let A,Ui, . . . ,U n be indeterminates over M, set U := (Ui, . . . , U n ) and 
denote by F : R ra+1 — > W 1 the nonlinear map defined by the right-hand side 
of ([7]). From the first n — 1 equations of (|7|) we easily see that, for a given 
positive value U\ = u%, the (positive) values of U 2 , ■ ■ ■ , U n ,A are uniquely 
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determined. Therefore, letting U\ vary, we may consider U2, ■ ■ ■ ,U n , A as 
functions of U\, which are indeed recursively defined as follows: 

U\{u\) := Ui, 
U 2 {ui) := ui + ^-gi(ui), 
U k+1 ( Ul ) := 2?7 A; (n 1 )-C/ fc _ 1 (n 1 ) + /i 2 5l ([/ fe (n 1 )), (2<fc<n-l), 

A( Ul ) := (\(U n -U n ^ l )(u l ) + ^g l {U n {u l )))/g 2 {U n (u 1 )). 

(8) 

Arguing recursively, one deduces the following lemma (cf. (2^, Remark 20]). 
Lemma 1. For any U\ > 0, the following assertions hold: 

i. (u k - c/fc-iXm) = h 2 (\ gi { Ul ) + E-=2 5i(^'K))) > 0, 

it. U k (u 1 )= Ul + h 2 ^g 1 (u 1 ) + Z k J Z 1 2(k-j)9i(U j (u 1 ))^ > 0, 

in- K - U'^im) = hWg'^ux) + Y,)Z\ g'^U^U^ux)) > 0, 
iv. U' k ( Ul ) = l + h 2 (^g[(u 1 ) + ZjZUk-j)9m(u 1 ))U' J (u 1 )) > 1, 

for 2 < k < n. 

For the proof of the next lema we need the following technical result 
Remark 2. Let fi,f 2 ,f3 £ C 2 (M>o) be positive functions such that 

• fi(x) > 0, 

• f 2 {x) > 0, &{x) > 0, 

• fz(x) > h{x), 

for all x > 0. Let F : IR>o — > R 6e i/ie function defined by 

FM >= /i(/ 2 (x))-/i(/ 3 (x)) 
Mx)-/ 8 (x) • 

T/ien > /or a// x > 0. 

Proof. Fix 2; > 0. By the definition of .F we have that 

F'(x)(f 2 (x)-f 3 (x)) = (f[W*))ffa) ~ fLW*MV>) 

-F(z)(f 3 (x) - &(z)) 
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holds. Prom the Mean Value Theorem, there exists £ € (/ 3 (x), fe(x)) with 
F(x) = /((£). Therefore 

F'(x)(/ 2 (x)-/ 3 (x)) = f/((/2(x))-/{(£))^(x)+f/i(£)-/{(/ 3 (x))V^(x). 



Since /{, / 2 and / 3 are strictly increasing functions, we conclude that F'(x) > 
0. □ 

Now we prove an important result for the existence and uniqueness of 
the solutions of (J7]) 

Lemma 3. For any U\ > 0, the following assertions hold: 

«■ (M)'^)<°> 

(ffrl)W>o, 

/or 2 < k < n. 

Proof. Let L, j : R>o — > R be the function defined by 



T i x gi(Uj) - gi(Ui) f x 



where 1 < i < j < n. From Remark [2] and Lemma [H we deduce that 

-giW-gxiUi 



L'iA»i) = r jj.-u- )^) >0 - (!,) 



By dSJ) we have 



——(mi) = 1 + — L 2 ,i(u 1 ). 
9i(Ui) 2 

Combining these identities with ([9]) we obtain ©, (|mj) and (|u>]) for 
k = 2. Now, arguing inductively, suppose that our statement is true for a 
given k > 2. From (J8j) we have: 

U k+ i-U k \ ( U k -U x \-i 

1+ i^^r + iw7r) ) {ui) - 
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Applying the inductive hypotheses we deduce that flmj) holds for k + 1. On 
the other hand, by Lemma [1] and ([8]) we have that 

51 



fe+i, 



L ("l)=(#^ («l) 



— +/t + £fc+l,fc (Ml), 



t/fe+i-C/i. . / c/i(f/ fc ) U k+ i-U) 



n? \(Vi)=[y T -rr + L k+ i :k — — ) [u x ) 

gi[U k+1 ) \U k+1 -Ui U k+1 -Ui' 



-1 



U k -U k - X , f/ fc -C/i 2 \-i C/fc+i-f/feV/ ^ 



ffi(^fc) * ' U k+ x-Ux 

gi(U k+ 1 ) f _ ^_( T U k+1 -Ui 
9i(Ui 



v 5i(^l) y 



^(w(^+£(*-^))M + i. 



hold. Combining the inductive hypothesis with ([9]) and (|mj) . for fe + 1, we 
conclude that © , (H3) and dwll hold for A; + 1 . □ 



2.2. Existence and uniqueness 

Let P : (]R>o) 2 -> K be the nonlinear map defined by 

P(a, Ul ) := \{U r ^ 1 {u 1 )-U n {u 1 ))-^g l {U n (u 1 )) + ag 2 {U n {u l )). (10) 

Observe that P(A, U\) = represents the minimal equation satisfied by the 
coordinates (a,u\) of any (complex) solution of the nonlinear system ([7]). 
Therefore, for fixed a > 0, the positive roots of P(a,Ui) are the values of 
ui we want to obtain. Furthermore, from the parametrizations (jHJ) of the 
coordinates u 2 , ■ ■ . , u n of a given solution (a, m, . . . , u n ) G (R>o) n+1 of (J7J) 
in terms of ui, we conclude that the number of positive roots of P(ct, U\) 
determines the number of positive solutions of ([7]) for such a value of a. 

Therefore, we analyze the existence of positive roots of the function 
P(a, U\) for values a > 0. Let g : M>o — > M>o be the function defined by 

g(x) := *(x). (11) 

92 



By Lemma [TJj3) we have that 
P{a,ui) = ag 2 {U n (ui)) - h(\gx(ux) + YJjZl 9\{Uj{ui)) + \gi{U n {u{))) 
> ag 2 {U n (ui)) - gi(U n (ui)) = g 2 (U n (ui))(a - g(U n (ui))) 
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holds for any u\ > 0. 

Suppose that g is surjective. Therefore, there exist it*, uf > such that 
g{U n {u\)) = a and g(U n (u**)) = 2a/h hold. From this choice of u\ and u*[* 
and the inequality above we deduce 

P(a,ut) > g 2 {U n {uX))(a-g{U n {u\))) = 0, 

P(a,nf) = \{U n ^{uT)-U n {ur))+g 2 (U n (ur))[a-^g{U n (uT))) 
= j^(l7 n _i(uj*)-l7 n («f))<0. 

Since P(A,U\) is a continuous function in (M>o) 2 , from the previous con- 
siderations we obtain the following result. 

Proposition 4. Fix a > and n € N. // £/ie function g of is surjective, 
then ([7]) has a positive solution with A = a. 

In order to establish the uniqueness, we prove that the homotopy path 
that we obtain by moving the parameter a in M > q is smooth. For this 
purpose, we show that the rational function A(U\) implicitly defined by the 
equation P(A, U\) = is decreasing. We observe that an explicit expression 
for this function in terms of U\ is obtained in ([8]). 

Theorem 5. Let A(U\) be the rational function of |2p. If the function g 
°f {HP is decreasing, then the condition A'(u\) < is satisfied for every 
ui € K >0 . 

Proof. Let U\, U2, ■ ■ ■ , U n ,A be the functions defined in ©. Observe that 
A can be rewritten as follows: 

A (TT \( Un ~ Un - 1 , k \ 

Taking derivatives with respect to U\, we have 

A' 'ITT \tt> ( Un ~ Un-l , h\ ( Un - U n -l \' 

Fix iti > 0. From Lemma[T]we see that g{U n ){u\) is positive. Furthermore, 
by Lemma El we have 

fU n -U n -i\' . 



9 



These remarks show that 



A> M < ( s ^„ W (^^ + f)) w _ ( s < (£/ „ W _A_) (ai) . (12) 

From Lemma [H we deduce that U' n (u\) and A{u\) are positive. Combining 
this affirmation with the monotonicity of g we deduce the statement of the 
theorem. □ 

Now we state and prove the main result of this section. 

Theorem 6. Let be given a > and n £ N. If the function g of ill]) is 
surjective and decreasing, then (jBJ) has a unique positive solution. 

Proof. Proposition U] shows that (jHJ) has solutions in (lR > o) n for any a > 
and any n € N. Therefore, there remains to show the uniqueness assertion. 

By Theorem[5l the condition A'{u\) < holds for every u± £ M>o- Argu- 
ing by contradiction, assume that there exist two distinct positive solutions 
. . . ,u n ), (ui, . . . ,u n ) £ (M>o) n of © for a. This implies that u% ^ U\ 
and A{u\) = A(ui), where A(Ui) is defined in ([8]). But this contradicts the 
fact that A'(ui) < holds in M>o, showing thus the theorem. □ 



3. Bounds for the positive solution 

In this section we obtain bounds for the positive solution of (|7|). More 
precisely, we find an interval containing the positive solution of ([7]) whose 
endpoints only depend on a. These bounds will allow us to establish an 
efficient procedure of approximation of this solution. 

Lemma 7. Let (a,u) G (M>o) n+1 be a solution of (|7|) for A = a. Then 

ug2{u n ) < g\(u n ). 

Proof. From the last equation of ([7]) for A = a and Lemma [DJI]) , we obtain 
the identity 

ag2{u n ) = h^-g^m) + gi(u 2 ) H \- gi(u n -i) + 

From the previous identity and Lemma [TJE]) we immediately deduce the 
statement of the lemma. □ 

From Lemma [7] we obtain the following corollary. 
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Corollary 8. Let (a,u) G (M >0 ) n+1 be a solution of flTJ for A = a. If the 

function g of ill]) is surjective and strictly decreasing, then 

Un < 5 _1 (a)- 

Let (a, u) G (R>o) n+1 be a solution of ([7]) for A = a. In the following 
lemma we obtain an upper bound of u n in terms of u± and a. 

Lemma 9. Let (a,u) G (K>o) ra+1 be a solution of ([7]) for A = a. If the 
function g of 177]) is surjective and strictly decreasing, then u n < e M u\ 
holds, with M := g'i{g~ l (a)) . 

Proof. Let (a, u) E (lR>o) ra+1 be a solution of ([7]) for A = a. Combining 
Lemma[JJJi|) and the Mean Value Theorem, we obtain the following identities 

Uk+i = Uk + h 2 (s±&l + gi(u2) H \-gi(u k )\ 

= u k + h 2 + g[(&)u 2 + ■■■+ g[^ k )u k ) 

for 1 < k < (n — 1), where £j G [0, nj] for 1 < i < k. Since g[ is an increasing 
function in M>o, combining Lemma [TJJi]) and Corollary [8l we obtain 

u k+ i < u k + h?(g' 1 (ui)ui H hsi^fc)^) 

< (l + V 1 (p- 1 (a))K = (l + fcM)« Jfc 

for 1 < A; < (n — 1). Arguing recursively, we deduce that 

u n < (l + Mh) n - l Ul < e M u x . 

This completes the proof. □ 

In our next lemma we obtain a lower bound of u\ in terms of a. 

Lemma 10. Let (a,u) G (M>o) ?1+1 be a solution of ([7J for A = a. If the 

function g of lill\) is surjective and strictly decreasing, then 

u\ > g^ 1 (aC(a)) 
holds, where C(a) > 1 is a constant such that 

lim C(a) = 1. 

a— >+oo 

Proof. Prom Lemma and Lemma [TJJi]) we deduce the inequalities 
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• 92{Un) < 52(e M Ui), 



• 91 (ui) < h(±gi (ui) + 5i(u 2 ) H h fll(«n-l) + 35i(m„)J = ag 2 (u n ). 

with M := g^(g -1 (a)). Combining both inequalities we obtain 

01 («l) < ag 2 {e M u 1 ) = a 92 ^ U ^ g 2 (u{). (13) 

92 (mi) 

Since (72 is an analytic function in x = and g 2 (x) 7^ for every x > 0, 
exists > 1 such that 

lim = e fcM . 

x-H)+ 02 (z) 

Combining this with Corollary [8l we deduce that there exists C(a) > such 
that 

1 < i^^L < C (a). 
32 («l) 

Furthermore, we can choose C(a) with 

lim C(a) = 1. 

a— >+oo 

Combining this remark with (|13p we obtain 

5i(ni) < aC(a)0 2 (Mi), 
which immediately implies the statement of the lemma. □ 

4. Numerical conditioning 

Let be given n £ N and a* > 0. In order to compute the positive solution 
of (J7]) for this value of n and A = a*, we shall consider (J7]) as a family of 
systems parametrized by the values a of A, following the positive real path 
determined by ([7j) when A runs through a suitable interval whose endpoints 
are a* and a*, where a* is a positive constant independent of h to be fixed 
in Section O 

A critical measure for the complexity of this procedure is the condition 
number of the path considered, which is essentially determined by the inverse 
of the Jacobian matrix of ([7]) with respect to the variables U\, . . . , U n , and 
the gradient vector of ([7]) with respect to the variable A on the path. In this 
section we prove the invertibility of such a Jacobian matrix, and obtain an 
explicit form of its inverse. Then we obtain an upper bound on the condition 
number of the path under consideration. 
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4-1. The Jacobian matrix 

Let F := F(A, U) : M n+1 -> W 1 be the nonlinear map denned by the 
right-hand side of ([7]). In this section we analyze the invertibility of the 
Jacobian matrix of F with respect to the variables U, namely, 



AAU) : =^(A,U): 



( 



\ 



■l 



-i r„ 



/ 



with Tx := 1 + \h 2 g[{Ui), r» := 2 + h 2 ^{Ui) for 2 < i < n - 1 and 

We start relating the nonsingularity of the Jacobian matrix J (a, u) with 
that of the corresponding point in the path determined by (|7|). Let (a, u) G 
1 be a solution of ([7]) for A = a. Taking derivatives with respect 
and substituting u\ for U\ we obtain the following tridiagonal 



to Ux in ((S 

system: 

/ ri(«i) 
-i 

V 



\ / 1 N 

t^(«i) 



1 



Tn(ux) J 



V hg 2 {U n (ux))A'{ux) ) 



For 1 < k < n — 1, we denote by A/% := Afc(^4, [/) the fcth principal minor of 
the matrix J (A, U), that is, the (k x fc)-matrix formed by the first k rows and 
the first k columns of J (A, U). By the Cramer rule we deduce the identities: 



hg 2 (U n (ux))A'(ux) 
det (J(a, u))U' k {ux) 



det (J (a, u)), 

hg 2 {U n (ux))A' (m) det (A fc _i(a, n)), 



(14) 
(15) 



for 2 < k < n. Suppose that a > and that the function g of (|lip is 
decreasing. Then Theorem [5] asserts that A'{u\) < holds. Combining this 
inequality with f)14[) we conclude that det (J(a, u)) < holds. Furthermore, 
by (fT5|) . we have 



L^(«!) = det(A fc _i(a,u)) (2 < Jfe < n). 



(16) 



Combining Remark [T lfn^ and (|16p it follows that det (A^a, u)) > holds 
for 1 < k < n — 1. As a consequence, we have that all the principal minors of 
the symmetric matrix A n _i(a,ii) are positive. Then the Sylvester criterion 
shows that A n _x(a,"u) is positive definite. These remarks allow us to prove 
the following result. 
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Theorem 11. Let (a,v) G (M >0 ) n+1 be a solution of for A = a. If the 

function g of ill]) is decreasing, then the matrix J(a, u) is invertible with 
det (J (a, u)) < 0. Furthermore, their (n-l)th principal minor is symmetric 
and positive definite. 

Having shown the invertibility of the matrix J(a, u) for every solution 
(a, u) £ (M>o) n+1 of (0, the next step is to obtain an explicit expression for 
the corresponding inverse matrices J -1 (a, u). For this purpose, we establish 
a result on the structure of the matrix J (a,u). 

Proposition 12. Let (a,u) 6 (M>o) n+1 be a solution of ([?]). If the function 
g of ill]) is decreasing, then the following matrix factorization holds: 



J (a, u) = 



where d{J) := det (J(a,u)) and u' k := U' k (ui) for 2 < k < n. 

Proof. Since J(a, u) is symmetric, invertible, tridiagonal and their (n — l)th 
principal minor is positive definite, the proof follows by a similar argument 
to that of [Ha, Proposition 25]. □ 

4-2. Upper bounds on the condition number 

From the explicit expression of the inverse of the Jacobian matrix J(A, U) 
on the points of the real path determined by (|7|), we can finally obtain esti- 
mates on the condition number of such a path. 

Let a* > and a* > be given constants independent of h. Suppose 
that the function g of (|lip is surjective and decreasing. Then Theorem [6] 
proves that (|7|) has a unique positive solution with A = a for every a in the 
real interval 1 := X(a*, a*) whose endpoints are a* and a*, which we denote 
by (ui(a), U2(ui(a)), . . . , U n {u\{a))). We bound the condition number 

K := max-fll^a^Hoo : a G X}, 

associated to the function ip : 1 — > M. n , (p(a) := (ui(a),U2(ui(a)), . . . , 
U n (ui(a))). 
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/ 1 



1 "2 
1 7? 





«2 


U n-1 


u 'n 


< 


< 


1 






dfj) 


dfj) ■ 





"3 
d(J) 



For this purpose, from the Implicit Function Theorem we have 

idF 



(—{a,<p(a)) 



dA 



(a,<p(a)) 



dF 



J (a,v?(a)) — (a, 99(a)) 



l^(a)Hc 



We observe that (dF/dA)(a,<p(a)) = ( 0, . . . , 0, -hg 2 (U n (ui(a))) ) holds. 
From Proposition [T2] we obtain 



l^(«)llc 



hg 2 (U n (u 1 (a))) / , , n< 

l,U 2 (ui{a)),... ,U n {ui{a)) 



det (J(a, tp(a))) 
Combining this identity with f|14|) . we conclude that 

1 



l^(«)|| 



1, U 2 (ui (a)),. ..,U' n (ui (a)) 



A'(ui(a)) 

From Lemma [H we deduce the following proposition. 

Proposition 13. Let a* > and a* > be given constants independent of 
h. Suppose that the function g of Ml}) is surjective and decreasing. Then 



I^HIIc 



K( Ul (a)) 



\A'( Ul (a))\ 
holds for a El. 

Combining Proposition [13] and (|12p we conclude that 

g (t7n(ui(q))) 



l^(«)Hc 



< 



a| 9 '(a„(«i (a)))|' 

Applying Lemma [TUl and Corollary [8] we deduce the following result. 

Theorem 14. Let a* > and a* > be given constants independent of h. 
Suppose that the function g of ill]) is surjective and g'{x) < holds for all 
x £ M>o. Then there exists a constant Ki(a*,a*) > independent of h with 



k < Ki(a*, a*). 
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5. An efficient numerical algorithm 

As a consequence of the well conditioning of the positive solutions of 
([7]), we shall exhibit an algorithm computing the positive solution of ([7]) for 



A = a* . This algorithm is a homotopy continuation method (see, e.g., 31 
§10.4], (29I . §14.3]) having a cost which is linear in n. 

There are two different approaches to estimate the cost of our procedure: 



using Kantorovich-type estimates as in [3l|, §10.4], and using Smale-type 



estimates as in 



29 . §14.3]. We shall use the former, since we are able to 
control the condition number in suitable neighborhoods of the real paths 
determined by ([7]). Furthermore, the latter does not provide significantly 
better estimates. 

Let a* > be a constant independent of h. Suppose that the following 
conditions hold: 

• g is surjective, 

• g'(x) < holds for all x > 0, 

• g"{x) > holds for all x > 0, 

where g is the function of (fTTj) . Then the path defined by the positive 
solutions of dTJ) with a € [a*, a*] is smooth, and the estimate of Theorem [T4l 
hold. Assume that we are given a suitable approximation of the positive 
solution ip(a*) of dTJ) for A = a*. In this section we exhibit an algorithm 
which, on input u^°\ computes an approximation of (p(a*). We recall that 
if denotes the function which maps each a > to the positive solution of 
([7J) for A = a. 

From Corollary [8] and Lemma \W\ we have that the coordinates of the 
positive solution of (jTJ tend to zero when a tends to infinity. Therefore, for 
a large enough, we obtain a suitable approximation of the positive solution 
([TJ for A = a*, and we track the positive real path determined by flTJ) until 
A = a*. 

In order to deal with a bounded set, we consider the change of variables 
B := 1/A. Then system (jJJ) for A = a* can be rewritten in terms of B as 
follows: 

' = -(U 2 -U 1 ) + ^g 1 (U 1 ), 

= -(U k+1 -2U k + U k . 1 ) + h 2 g 1 (U k ), (2<fc<n-l), 
= -{U n . l -U n )-hB- l g 2 (U n ) + ^g l (U n ) 

(17) 
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for B = /?*, where /?* := 1/a*. 

Let < /3* < (3* be a constant independent of /i to be determined. Fix 
P E By Corollary [8] it follows that (f(a) is an interior point of the 

compact set 

Kp ■= Hoc <2g- 1 (l/P)}, 

where cp : [/?*,/?*] — > R n is the function which maps each (3 E [/?*,/3*] to the 
positive solution of (|17j) for B = (3, namely 

¥>(/3) := (ttiGS), . . . ,u n ((3)) := U 2 {u x {fi% U n ( Ul (f3)))- 

First we prove that the Jacobian matrix Jfi(u) := (dF /dU)(f3,u) is in- 
vertible in a suitable subset of Kr. Let u E M ri and v E M ri be points 
with 

\\u - vKpOIU < 5/3, lb - ¥>(/3)||oo < S/3, 

where 5p > is a constant to be determined. Note that if 5p < 

then u E If a and t> E if^. By the Mean Value Theorem, we see that the 

entries of the diagonal matrix Jp{u) — Jp(v) satisfy the estimates 

(Jp(u)- Jp(v)) u \ < 2h 2 g'{(2g- 1 (l/P))5 p (l<i<n-l), 

{Jp{u)-J p {v)) nn \ < 2h^{&{2g-\l/f3))/p,gi{(2g-\l/f3))}tp- 

By Theorem [11] and Proposition [12] we have that the matrix Jiptp) '■= 
Jp(<p(/3)) = (dF/dU)(P,<p(p)) is invertible and 

M ^ gMg)g(t^)) uKumwMP)) 

m j fc=m tr {M} + det(^ (/3) ) 

holds for 1 < i, j < n. According to Lemma[TJ we have U' n {ui{l3)) > • • • > 
U^iui (/3)) > 1. These remarks show that 



< 

oo 



J tp(p)( J p( u ) - J ^ v )) 

-^( 2 + ~ IdetfX^ ) ( 18 ) 
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where rfp := 2max{g'{(2g- 1 (1/ /3)), g'^(2g- 1 (l/ (3)) / /?}. Since B = l/A, from 
Theorem [5] we deduce that B'(ui) = - A' (ui) / A 2 (ui) > for all x > 0. 
Combining this assertion with ()14|) . we obtain the following identity: 



hu' n { Ul m _ u' n (um)B 2 (um) 

|det(J^))| B'( Ul (f3))g 2 (u n ((3)) 
From ()12p . we have that 



|det(J^ (/3) )| B'(ui{p))g 2 {u n (f3)) ~ \g'(u n (/3))\g 2 (u n ((3))' { ' 

Combining this inequality with the definition and the monotonicity of g, we 
deduce 

\g'{u n {f3))\g 2 {u n (P)) = giMPMMP)) ~ g / iK(/3))ff 2 K(/3)) 

g(un(P)) gi(u n (P)) 

Since g\ and g 2 are analytic functions in x = and 51(0) = g 2 {Q) = 0, there 
exists r > such that, for all |x| < r, we have that 

51 0*0 = ^2c k x k ,g 2 (x) = ^d k x k , 

k=p k=q 

with c p 7^ and d g 7^ 0, where p and (7 are positive integers greater than 1. 
Hence, the following identity holds: 

hmn _g' 1 (,x)g 2 ( x )\ =l _ P __ 



x->o \ gi(x)g' 2 (x) 

Taking into account that u n (/3) € (0,s _1 (l//3*)] holds for all /3 e [/3*,/3*], 
we conclude 

7 — t^tt > g 2 {u n {(3))(l - p , 

where is a constant which depends only on f3* . Combining the last in- 
equality with (fT8j) and (fT9|) . we obtain 



< 2^ 1 + - 



g>{u n ({3))\g 2 {u n {p)) 
B( Ul (f3)) \ 
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From Lemma PTOl and Corollary [HJ there exists a constant C(/3) > 1 such that 



^)(w-w)L<^(i+-^ 



g' 2 {g-\C{P)/pm-P*) 
Furthermore, the following condition holds: 

lim C(B) = 1. 
Finally, since g'{x) < holds for all x > 0, we have 

p _ cm 

g' 2 (g^(C(P)/P))(l-p*) 



(20) 



< 



g' 2 (g~HC(/3)//3))g(g-^C(f3)/f3))(l ~ P*) 
C(J3) 



A{srHC{P)IPm-fY 

Combining this inequality with (I20p . we deduce that 

2^(0* + l)C(/3) 



< 



rfy, (21) 



g[(g-HC((3)/mi-p*) * 
with 0* := (1 — p*)g' 1 (g~ 1 (l//3*)). Hence, defining 5p in the following way: 



mm 



g ;(g- 1 (C(/3)//3))(l-p*) 
8^(0* + l)C(/3) 



we obtain 



J v tf))[M u ) ~ M v ) 



< -. 

oo 4 



(22) 



(23) 



In particular, for v = (f(/3), this bound allows us to consider Jp(u) as a 
perturbation of J V (B)- More precisely, by a standard perturbation lemma 
(see, e.g., {3l|, Lemma 2.3.2]) we deduce that Jp(u) is invertible for every 
u G Bsp((p(f3)) fl -fTg and we obtain the following upper bound: 



Jp(u) 1 J ip (f3) 



< -. 

oo 3 



(24) 



In order to describe our method, we need a sufficient condition for the 
convergence of the standard Newton iteration associated to ([7]) for any f3 € 
[/?*,/?*]. Arguing as in 3jJ, 10.4.2] we deduce the following remark, which in 
particular implies that the Newton iteration under consideration converges. 
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Remark 15. Set 5 := min{£g : (3 € [/?*, /3*]}. Fix (3 G [/?*,/?*] and consider 
the Newton iteration 



u (k+i) = u (k) _ jpfaMylpfa „(*)) (k>0), 
starting at G JTg. // - ^(/3)||oo < S, then 



i [k) -m\\oc<^ 



holds for k > 0. 



Now we can describe our homotopy continuation method. Let (3q := 
/3* < Pi < ■ ■ ■ < [3n := f3* be a uniform partition of the interval [/3*,/3*], 
with N to be fixed. We define an iteration as follows: 

u (k+l) = u (k) _ jp k ( u (k))-i F (p k , U W) (0<k<N-l), (25) 

u (N+k+l) = u (N+k)_ Jf) .( u (N+k)yl F (p* iU (N+k)) (yfc > 0) _ (26) 

In order to see that the iteration (|25p - (j26p yields an approximation of 
the positive solution ip(f3*) of (fT7|l for B = /3*, it is necessary to obtain 
a condition assuring that (|25p yields an attraction point for the Newton 
iteration (I26D . This relies on a suitable choice for N, which we now discuss. 

By Theorem 1141 we have 

|U < maxdl^C^Hoo :j8e /T]} |/3 m 

for < i < JV — 1, where ki is an upper bound of the condition number 
independent of h. Thus, for N := \3/3*ki/S~\ + 1 = 0(1), by the previous 
estimate we obtain the following inequality: 

MA + i)-^(A)l!oo<^ (27) 

for < i < N — 1. Our next result shows that this implies the desired result. 

Lemma 16. Set N := \3f3*Ki/5~\ + 1. Then, for every u^ ' with \\u^ — 
< p(P*)\\od < $, the point v,( N ' defined in &25\) is an attraction point for the 
Newton iteration 
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Proof. By hypothesis, we have ||«^ — </?(/?*) ||oo < S. Arguing inductively, 
suppose that \\v\ k > — <p(Pk)\\oo < $ holds for a given < k < N. By 
Remark [15] we have that vS k ^ is an attraction point for the Newton iteration 
associated to (|17p for B = /3k- Furthermore, Remark 1151 also shows that 
|| n (fc+i) _ (ptfjWn < 6/3 holds. Then 

Hu^-^fc+Olloo < lk (fe+1) -^(A)||oc + ||^)-^+l)||oo 

< y + y < 5, 

where the inequality ||</?(/3fc+i) — <p(Pk)\\oc < S/3 follows by (j27|) . This com- 
pletes the inductive argument and shows in particular that u^ N ' is an at- 
traction point for the Newton iteration (|26p . □ 



Next we consider the convergence of (|26p . starting with a point u^ N ' sat- 
isfying the condition \\u^ N > — ip((3*)\\oo < 8 < S/3*- Combining this inequality 
with (|22p we deduce that € K a * . Furthermore, we see that 

||«(^+ 1 )-^( ) 8*)|| oo = ||uW-J^(«W)- 1 F( j 9*,uW)- V C8*)|| 0o 

= J p , (uW)- 1 ^. (#)) («W -<ptf*)) -F(0*, 08*, ))) 

OO 

|^)(V(^ W )(« W -^(/3*))-^(/3*^ (JV) )+^^(/3*))) 

(28) 

where £ is a point in the segment joining the points vP^ and (p(/3*). Com- 
bining (f2Tj) and we deduce that 



OO 
OO 1 



l (^+i)_^* ) || oo < g||j^. ) (j /J .(tiW)-j /J .(0)|| oo ^ 



holds, with c:= (2^.(0* + l)C(f3*))/(g[(g- l (C((3*)/f3*))a ~ P*))- By an 
inductive argument we conclude that the iteration (I26p is well-defined and 
converges to the positive solution (p(f3*) of ()1T|) for B = f3*. Furthermore, 
we conclude that the point u^ N+k \ obtained from the point u^ N ' above after 
k steps of the iteration ([26]) . satisfies the estimate 



||n^)-^)||oo<c(fV) 2fc <<^) 2 '> 

with c := 3/4c. Therefore, in order to obtain an e-approximation of <£>(/?*), 
we have to perform log 2 log 3 (3/4ce) steps of the iteration ([26]) . Summarizing, 
we have the following result. 
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Lemma 17. Let e > be given. Then, for every G (R>o) n satisfying 
the condition \\u^ N ' — (p((3*)\\ OD < 5, the iteration 126}) is well-defined and 
the estimate \\u^ N+k ^ - v?(/3*)||oo < £ holds for k > log 2 log 3 (3/4ce). 

Let £ > 0. Assume that we are given G (M>o) n such that \\u^ — 
¥>(/3*)||oo < 5 holds. In order to compute an e-approximation of the positive 
solution (f((3*) of (|17p for B = (3*,we perform N iterations of ()25[) and k$ := 
[log 2 log 3 (3/4ce)] iterations of ([26]) . From Lemmas [TBI and [T71 we conclude 
that the output u^ N+k °' > of this procedure satisfies the condition — 
y ) ||oo < £• Observe that the Jacobian matrix Jp{u) is tridiagonal for 
every (3 G [/3*,/3*] and every u G Kp. Therefore, the solution of a linear 
system with matrix Jp{u) can be obtained with 0(n) flops. This implies 
that each iteration of both (|25p and (|26p requires O(n) flops. In conclusion, 
we have the following result. 

Proposition 18. Let e > and G (M >0 ) n loii/i ||ii (0) - <p(P*)\\oo < 
5 be given, where 5 is defined as in Remark 1151 Then the output of the 
iteration i25\) -[2S \) is an e-approximation of the positive solution (p(ft*) of 
pl\ ) for B = (3* . This iteration can be computed with 0{Nn + k$n) = 
0{n log 2 log 2 (l/e)) flops. 

Finally, we exhibit a starting point G (M>o) n satisfying the condi- 
tion of Proposition [181 Let /3* > be a constant independent of h to be 
determined. We study the constant 

6:=wm{5f,:P€\p m ,P*]}, 

where 

^- mm i 8^(0* + l)C(/3) ' 5 

Since gi and g 2 are analytic functions in x = 0, in a neighborhood of G K n , 
we can rewrite ^ as follows: 

r,p = 2m a x{g'{(2g- 1 (l/t3)),g^2g- 1 (l/(3))g(g- 1 (l/m 
= 2max{5 1 (/3),S 2 (/3)}( 5 - 1 (l//3)) p - 2 , 

where p is the multiplicity of as a root of g\ and Si is an analytic function 
in x = such that limg-^o 7^ for i = 1,2. Taking into account that 

j3 G (0, j3*\ holds, we conclude that there exists a constant rf > 0, which 
depends only on /?*, with 

^<2r ? *( 5 - 1 (l//3)r 2 
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for all /3 G (0, (3*\. Moreover, with a similar argument we deduce that there 
exists a constant •&* > 0, which depends only on B*, such that 



5/3 > min 
We claim that 



f **(W) ( g^jCWm y- 1 Aa-UUB) (29) 



16t7*(6>* + 1)C(B) \ g- l {l/P) 

i im i&wm . ,-. ,30) 

/?->o+ 5 _1 (l//3) 

In fact, since we have C(/3) > 1 and g~ l is decreasing, it follows that 

9-HC(B)/B) 

On the other hand, there exists £ G (l/(3,C(B)/B) with 

9-HC(B)/B) ff- 1 (l//?) + (g- 1 y(Q((C(/3)-l)//3) 

Since g\ and 52 are analytic functions in x = and lirng_> + = 1, we 
see that 

hm 9{g-\l/P)) (c(g) _ 1) _ Q 



Combining (j31j) . (|32|) and this inequality we inmediately deduce (|30|) . 

Combining ([29]) with ([30j) we conclude that there exists a constant C* G 
(0, 1], which depends only on B*, with 

Sp > C*g-\l/B). 

Therefore, 

5 = wm{8p : /3 G [&,/?*]} > C*g-\l/B*). (33) 
From Corollary [8] and Lemma [TC)| we have 

<P(B*) elg-HctfJ/B^g-^l/B.)]". 

Furthermore, by (f30j) . we deduce that 

g-^CtfJ/B.) 



-)^ 1 (l//3*)<CV 1 (V/3*) (34) 
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holds for /3* > small enough. Combining this with (|33p . we conclude that 

||n - ^(/3*)||oc < - g~\C{^)/^) < 5 (35) 

holds for all u G [cT 1 ^/ 3 *)// 3 *), ^H 1 // 9 *)]™- T1 ™s, let /3* < 0* satisfy 
dMD - Then, for any in the hypercube Ig' 1 (C(P*)/ P*) , g' 1 (1/ P*)] n , th e 
inequality 

||« (0) -^*)||oo <5 
holds. Therefore, applying Proposition [TBI we obtain our main result. 

Theorem 19. Lete>0 be given. Then we can compute an e- approximation 
of the positive solution of (T7\ ) for B = j3* with 0(n log 2 log 2 (l/e)) flops. 
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